clear
input str12 region policy_year
"广东"    2012
"吉林"    2013
"北京"    2013
"浙江"    2013
"黑龙江"  2013
"江苏"    2013
"安徽"    2013
"山东"    2013
"湖南"    2013
"四川"    2013
"甘肃"    2013
"天津"    2014
"上海"    2014
"内蒙古"  2014
"河南"    2014
"湖北"    2014
"贵州"    2014
"云南"    2014
"重庆"    2014
"河北"    2014
"宁夏"    2014
"陕西"    2015
"山西"    2015
"广西"    2015
end

save "policy_years.dta", replace

// For policies published in November or December of that year, the implementation year is set to the following year

use "data.dta" ,clear       // *************************use data.xlsx*************************************

* Rename variables to English abbreviations
rename 地区 region
rename 年份 year
rename 基金会单位数个 found_num
rename 社会团体单位数个 social_org_num
rename 民办非企业单位数个 priv_npe_num
rename 社会组织总数 total_social_org
rename 社区居委会村委会合计单位数个自治组织 total_comm_num
rename 民政部门直接接收捐赠情况捐赠款数额万元 direct_donation            
rename 民政管理事务民间组织管理万元 ngo_manage_exp
rename 地区生产总值亿元 rgdp
rename 劳动报酬占GDP比重 labor_gdp_ratio
rename 全体居民人均消费支出 cons_exp_all
rename 城镇化率 urban_rate
rename 人均GDP增长率 gdp_growth
rename 私营企业就业人数万人 private_emp_total
rename 民办非企业年末职工人数人            priv_npe_emp
rename 基金会年末职工人数人                found_emp
rename 社会团体年末职工人数人              social_org_emp
rename 总人口人                            pop_total_persons
rename 社区服务中心年末职工人数人          comm_srv_emp
rename 社会捐赠款物合计亿元        social_don_total_yi
rename 各类社会组织社会捐赠款亿元  so_don_cash_yi
rename F社区服务中心年末职工人数人         f_comm_srv_emp
rename 年末常住人口万人 pop_total

foreach v in ///
    priv_npe_emp found_emp social_org_emp ///
    pop_total_persons f_comm_srv_emp comm_srv_emp {

    capture confirm string variable `v'
    if !_rc {
        replace `v' = "" if trim(`v')==""        // Replace whitespace-only strings with empty strings
        destring `v', replace ignore(", ")       // Convert string to numeric, ignoring commas and spaces
    }
}


foreach v in ///
    social_don_total_yi so_don_cash_yi {

    capture confirm string variable `v'
    if !_rc {
        replace `v' = ustrtrim(`v')                 // Trim Unicode whitespace from strings
        replace `v' = "" if `v'==""                 // Handle empty strings
        destring `v', replace ignore(", ， ")       // Convert to numeric, ignoring both English and Chinese commas
    }
}

isid region year, sort
list region year f_comm_srv_emp if missing(f_comm_srv_emp)
sort region year
by region: replace f_comm_srv_emp = (f_comm_srv_emp[_n-1] + f_comm_srv_emp[_n+1]) / 2 ///
    if missing(f_comm_srv_emp) ///
    & year[_n-1] == year-1 ///
    & year[_n+1] == year+1
count if missing(f_comm_srv_emp)
list region year f_comm_srv_emp if region=="西藏" & inrange(year, 2014, 2018), sep(0)

save "date_full.dta", replace

drop if year > 2016
save "date_2016.dta", replace

// Merge main dataset with policy mapping dataset
merge m:1 region using "policy_years.dta"
tab _merge
drop _merge

// Generate treatment variable
gen treatment = (year >= policy_year) if !missing(policy_year)
replace treatment = 0 if missing(treatment)
replace treatment = 1 if treatment >= 1
replace treatment = 0 if treatment < 1
label define treatlbl 0 "Control" 1 "Treated"
label values treatment treatlbl
sort region year

// Verification
tabulate treatment

gen private_emp_ratio = private_emp_total/pop_total
gen ln_direct_donation = ln(direct_donation)
gen ln_ngo_manage_exp = ln(ngo_manage_exp)
gen ln_private_emp_total   =ln(private_emp_total )
gen ln_cons_exp_all   =ln(cons_exp_all )
gen ln_pop = ln(pop_total_persons)
gen ln_srve = ln(f_comm_srv_emp)
encode region, generate(region_num)

* Overall descriptive statistics
estpost tabstat total_social_org found_num social_org_num priv_npe_num ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve ln_ngo_manage_exp social_don_total_yi so_don_cash_yi, statistics(N mean sd min p50 max) columns(statistics)

* Descriptive statistics by group
estpost tabstat total_social_org found_num social_org_num priv_npe_num ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve ln_ngo_manage_exp social_don_total_yi so_don_cash_yi, by(treated) statistics(N mean) columns(statistics)


// Part II: Staggered DID Regression Results - 4 Models

sort region year

reghdfe total_social_org treatment ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

 reghdfe found_num treatment ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

 reghdfe social_org_num treatment ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

reghdfe priv_npe_num treatment ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

// Robustness check
* 1) Calculate total (ensure found_emp has no missing values)
gen total_org_emp = priv_npe_emp + found_emp + social_org_emp

// Robustness check 2006-2016
// Total number of employees: total_org_emp
reghdfe total_org_emp treatment ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)
// Foundation employees at year-end (persons): found_emp
reghdfe found_emp treatment ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)
// Social organization employees at year-end (persons): social_org_emp
reghdfe social_org_emp treatment ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)
 // Private non-enterprise unit employees at year-end (persons): priv_npe_emp
reghdfe priv_npe_emp treatment ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

preserve
use "date_full.dta" ,clear
// Robustness check: the following four tables use 2006-2022 data, with NGO counts as the dependent variable
// private_emp_total (private enterprise employment in 10,000 persons) is only available until 2019
// labor_gdp_ratio (labor compensation as % of GDP) is only available until 2017
// ln_direct_donation (direct donations received by civil affairs departments in 10,000 yuan) is only available until 2017
// Therefore, these three control variables are not included in these four models
reghdfe total_social_org treatment total_comm_num ln_cons_exp_all urban_rate gdp_growth ln_pop ln_srve, absorb(year region)
reghdfe found_num treatment total_comm_num ln_cons_exp_all urban_rate gdp_growth ln_pop ln_srve, absorb(year region)
reghdfe social_org_num treatment total_comm_num ln_cons_exp_all urban_rate gdp_growth ln_pop ln_srve, absorb(year region)
reghdfe priv_npe_num treatment total_comm_num ln_cons_exp_all urban_rate gdp_growth ln_pop ln_srve, absorb(year region)
restore

 
preserve
use "date_full.dta" ,clear
// Robustness check: the following four tables use 2006-2022 data, with full-time employees as the dependent variable
// private_emp_total (private enterprise employment in 10,000 persons) is only available until 2019
// labor_gdp_ratio (labor compensation as % of GDP) is only available until 2017
// ln_direct_donation (direct donations received by civil affairs departments in 10,000 yuan) is only available until 2017
// Therefore, these three control variables are not included in these four models
reghdfe total_org_emp treatment   total_comm_num  ln_cons_exp_all urban_rate gdp_growth  ln_pop ln_srve, absorb(year region)
reghdfe found_emp treatment   total_comm_num  ln_cons_exp_all urban_rate gdp_growth  ln_pop ln_srve, absorb(year region)
reghdfe social_org_emp treatment   total_comm_num  ln_cons_exp_all urban_rate gdp_growth  ln_pop ln_srve, absorb(year region)
reghdfe priv_npe_emp treatment   total_comm_num  ln_cons_exp_all urban_rate gdp_growth  ln_pop ln_srve, absorb(year region)
restore





// Part III: Event Study Analysis
// Event study variables
// Generate treatment group indicator (if not already generated)
gen treated = !missing(policy_year)
label define treatlbl_new 0 "Control" 1 "Treated"
label values treated treatlbl_new

// Generate event time variable: years relative to policy implementation year
gen event_time = year - policy_year

// Define event time window, e.g., -5 to +5
gen event_time_window = event_time

// Limit event time to range -5 to +5
//replace event_time_window = -5 if event_time_window < -5
//replace event_time_window = 4 if event_time_window > 4

// Generate event time dummy variables, avoiding negative signs in variable names
forvalues i = -9/4 {
    // Generate different suffixes based on the value of i
    if `i' < 0 {
        local suffix = "m" + string(abs(`i'))
    }
    else if `i' > 0 {
        local suffix = "p" + string(`i')
    }
    else {
        local suffix = "0"
    }
    
    // Generate dummy variable
    gen event_time_`suffix' = (event_time_window == `i') & treated
}

// Check generated variables
describe event_time_*

* Plot combined event study graph	
* Set unified plotting parameters
local graph_opts ///
    vertical ///
    yline(0, lcolor(gs8) lpattern(dash)) ///
    graphregion(color(white) margin(small)) ///
    plotregion(margin(tiny)) ///
    omitted base ///
    xtitle("Event Time", size(small)) ///
    legend(cols(2) size(small) region(lcolor(white)))
	
* Run four regressions (excluding event_time_m2 as reference period)
* Regression 1: total_social_org
reghdfe total_social_org event_time_m9 event_time_m8 event_time_m7 event_time_m6 event_time_m5 event_time_m4 event_time_m3 event_time_m1 event_time_0 event_time_p1 event_time_p2 event_time_p3 event_time_p4 ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)
estimates store total

* Regression 2: found_num
reghdfe found_num event_time_m9 event_time_m8 event_time_m7 event_time_m6 event_time_m5 event_time_m4 event_time_m3 event_time_m1 event_time_0 event_time_p1 event_time_p2 event_time_p3 event_time_p4 ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)
estimates store found

* Regression 3: social_org_num
reghdfe social_org_num event_time_m9 event_time_m8 event_time_m7 event_time_m6 event_time_m5 event_time_m4 event_time_m3 event_time_m1 event_time_0 event_time_p1 event_time_p2 event_time_p3 event_time_p4 ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)
estimates store social

* Regression 4: priv_npe_num
reghdfe priv_npe_num event_time_m9 event_time_m8 event_time_m7 event_time_m6 event_time_m5 event_time_m4 event_time_m3 event_time_m1 event_time_0 event_time_p1 event_time_p2 event_time_p3 event_time_p4 ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)
estimates store private

* Plot combined graph: Total NGOs (NGO_total), Foundation count (FOUND_num), Social organization count (SA_num), and Private non-enterprise unit count (PNU_num)
coefplot ///
    (total, mcolor(black) msymbol(circle) ciopts(recast(rcap) lcolor(black) lwidth(thin)) offset(-0.2)) ///
    (social, mcolor(black) msymbol(square) ciopts(recast(rcap) lcolor(black) lwidth(thin)) offset(-0.067)) ///
    (private, mcolor(black) msymbol(triangle) ciopts(recast(rcap) lcolor(black) lwidth(thin)) offset(0.067)) ///
    (found, mcolor(black) msymbol(diamond) ciopts(recast(rcap) lcolor(black) lwidth(thin) yaxis(2)) offset(0.2) yaxis(2) ci(95)), ///
    keep(event_time_*) ///
    coeflabels(event_time_m9="-9" event_time_m8="-8" event_time_m7="-7" event_time_m6="-6" ///
               event_time_m5="-5" event_time_m4="-4" event_time_m3="-3" ///
               event_time_m1="-1" event_time_0="0" event_time_p1="1" event_time_p2="2" ///
               event_time_p3="3" event_time_p4="4") ///
    ytitle("Treatment Effect (Total NGOs, Social orgs, Private NPOs)", size(small)) ///
    ytitle("Treatment Effect (Foundations)", size(small) axis(2)) ///
    ylabel(0(5000)25000, labsize(small) angle(0)) /// Left axis scale
    ylabel(-50(50)400, labsize(small) angle(0) axis(2)) /// Right axis scale
    legend(order(2 "Total NGOs" 4 "Social orgs" ///
                 6 "Private NPOs" 8 "Foundations")) ///
    levels(95) ///
    `graph_opts'	
	
// Part IV: Placebo Test

preserve
* First equation - total_social_org  
reghdfe total_social_org treatment ln_direct_donation  total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

permute treatment _b[treatment], reps(500) seed(12345)  saving("simulations1.dta", replace): ///
reghdfe total_social_org treatment ln_direct_donation  total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

use "simulations1.dta", clear
dpplot *pm*1, xline(2198.907, lpattern(dash) lcolor(black) extend) ///
    msize(small) mcolor(black) msymbol(circle_hollow) ///
    xtitle("Treatment Effect") ytitle("Density") ///
    title("Total NGOs", size(medium)) ///
    scheme(s2mono) ///
    ylabel(,nogrid) xlabel(,nogrid) ///
    graphregion(color(white)) bgcolor(white) ///
    xscale(range(-2500 2500)) ///
    saving("placebo1.gph", replace)
restore

preserve
* Second equation - found_num  
reghdfe found_num treatment ln_direct_donation  total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

permute treatment _b[treatment], reps(500) seed(12345) saving("simulations_found.dta", replace): ///
reghdfe found_num treatment ln_direct_donation  total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

use "simulations_found.dta", clear
dpplot *pm*1, xline(24.583, lpattern(dash) lcolor(black) extend) ///
    msize(small) mcolor(black) msymbol(circle_hollow) ///
    xtitle("Treatment Effect") ytitle("Density") ///
    title("Foundations", size(medium)) ///
    scheme(s2mono) ///
    ylabel(,nogrid) xlabel(,nogrid) ///
    graphregion(color(white)) bgcolor(white) ///
    xscale(range(-30 30)) ///
    saving("placebo2.gph", replace)
restore

preserve
* Third equation - social_org_num 
reghdfe social_org_num treatment ln_direct_donation  total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

permute treatment _b[treatment], reps(500) seed(12345) saving("simulations_social.dta", replace): ///
reghdfe social_org_num treatment ln_direct_donation  total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

use "simulations_social.dta", clear
dpplot *pm*1, xline(935.321, lpattern(dash) lcolor(black) extend) ///
    msize(small) mcolor(black) msymbol(circle_hollow) ///
    xtitle("Treatment Effect") ytitle("Density") ///
    title("Social orgs", size(medium)) ///
    scheme(s2mono) ///
    ylabel(,nogrid) xlabel(,nogrid) ///
    graphregion(color(white)) bgcolor(white) ///
    xscale(range(-1200 1200)) ///
    saving("placebo3.gph", replace)
restore

preserve
* Fourth equation - priv_npe_num  
reghdfe priv_npe_num treatment ln_direct_donation  total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

permute treatment _b[treatment], reps(500) seed(12345) saving("simulations_priv.dta", replace): ///
reghdfe priv_npe_num treatment ln_direct_donation  total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

use "simulations_priv.dta", clear
dpplot *pm*1, xline(1239.003, lpattern(dash) lcolor(black) extend) ///
    msize(small) mcolor(black) msymbol(circle_hollow) ///
    xtitle("Treatment Effect") ytitle("Density") ///
    title("Private NPOs", size(medium)) ///
    scheme(s2mono) ///
    ylabel(,nogrid) xlabel(,nogrid) ///
    graphregion(color(white)) bgcolor(white) ///
    xscale(range(-1300 1300)) ///
    saving("placebo4.gph", replace)
restore
		
* Combine the four saved graphs
graph combine placebo1.gph placebo2.gph placebo3.gph placebo4.gph, ///
    rows(2) cols(2) ///
    title("Placebo Tests for Treatment Effects", size(medium)) ///
    scheme(s2mono) ///
    graphregion(color(white)) ///
    note("Notes: Dashed lines indicate the actual treatment effects", size(small)) ///
    xsize(10) ysize(8) ///
    saving("combined_placebo.gph", replace)

	
// Bacon Decomposition

* Replace Chinese region names with English names
replace region = "Shanghai" if region == "上海"
replace region = "Yunnan" if region == "云南"
replace region = "Inner Mongolia" if region == "内蒙古"
replace region = "Beijing" if region == "北京"
replace region = "Jilin" if region == "吉林"
replace region = "Sichuan" if region == "四川"
replace region = "Tianjin" if region == "天津"
replace region = "Ningxia" if region == "宁夏"
replace region = "Anhui" if region == "安徽"
replace region = "Shandong" if region == "山东"
replace region = "Shanxi" if region == "山西"
replace region = "Guangdong" if region == "广东"
replace region = "Guangxi" if region == "广西"
replace region = "Xinjiang" if region == "新疆"
replace region = "Jiangsu" if region == "江苏"
replace region = "Jiangxi" if region == "江西"
replace region = "Hebei" if region == "河北"
replace region = "Henan" if region == "河南"
replace region = "Zhejiang" if region == "浙江"
replace region = "Hainan" if region == "海南"
replace region = "Hubei" if region == "湖北"
replace region = "Hunan" if region == "湖南"
replace region = "Gansu" if region == "甘肃"
replace region = "Fujian" if region == "福建"
replace region = "Tibet" if region == "西藏"
replace region = "Guizhou" if region == "贵州"
replace region = "Liaoning" if region == "辽宁"
replace region = "Chongqing" if region == "重庆"
replace region = "Shaanxi" if region == "陕西"
replace region = "Qinghai" if region == "青海"
replace region = "Heilongjiang" if region == "黑龙江"

* Check replacement results
tab region, nol

panelview total_social_org treatment ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, ///
    i(region) t(year) type(treat) title("Treatment status") ytitle("Provincial administrative regions") xtitle("Year") name(cohort, replace) bytiming mycolor(Greys)

ddtiming total_social_org treatment,i(region_num) t(year)


// Interaction Effects
// Variable centering
egen mean_ln_ngo = mean(ln_ngo_manage_exp)  // Calculate mean of ln_ngo_manage_exp
egen mean_social_don_total_yi = mean(social_don_total_yi)  // Calculate mean of social_don_total_yi
egen mean_so_don_cash_yi = mean(so_don_cash_yi)  // Calculate mean of so_don_cash_yi

generate c_ln_ngo = ln_ngo_manage_exp - mean_ln_ngo    // Center ln_ngo_manage_exp
generate c_social_don_total_yi = social_don_total_yi - mean_social_don_total_yi  // Center social_don_total_yi
generate c_so_don_cash_yi = so_don_cash_yi - mean_so_don_cash_yi  // Center so_don_cash_yi

* Model 1: Baseline model (no interaction terms)
reghdfe total_social_org treatment ln_direct_donation total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

* Model 2: Pairwise interaction terms 
* Treatment × c_so_don_cash_yi interaction  c_social_don_total_yi

reghdfe total_social_org treatment so_don_cash_yi ln_direct_donation  total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

reghdfe total_social_org c.c_treatment##c.c_so_don_cash_yi  ln_direct_donation total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

reghdfe found_num        c.c_treatment##c.c_so_don_cash_yi  ln_direct_donation total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

reghdfe social_org_num   c.c_treatment##c.c_so_don_cash_yi  ln_direct_donation total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

reghdfe priv_npe_num      c.c_treatment##c.c_so_don_cash_yi  ln_direct_donation total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)


reghdfe total_social_org  treatment social_don_total_yi ln_direct_donation  total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

reghdfe total_social_org c.c_treatment##c.c_social_don_total_yi  ln_direct_donation total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

reghdfe found_num        c.c_treatment##c.c_social_don_total_yi  ln_direct_donation total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

reghdfe social_org_num   c.c_treatment##c.c_social_don_total_yi  ln_direct_donation total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

reghdfe priv_npe_num      c.c_treatment##c.c_social_don_total_yi  ln_direct_donation total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)



* Treatment × ln_ngo_manage_exp interaction

reghdfe total_social_org treatment ln_ngo_manage_exp ln_direct total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)


reghdfe total_social_org c.c_treatment##c.c_ln_ngo c_ln_direct total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

reghdfe found_num c.c_treatment##c.c_ln_ngo c_ln_direct total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

reghdfe social_org_num c.c_treatment##c.c_ln_ngo c_ln_direct total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)

reghdfe priv_npe_num c.c_treatment##c.c_ln_ngo c_ln_direct total_comm_num  labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve, absorb(year region)





	
	


	
// Part VI: Quantile Staggered DID
// First encode the region variable  

foreach q in 10 20 30 40 50 60 70 80 90 {
    qreg2 found_num treatment ln_direct_donation  ///
         total_comm_num  labor_gdp_ratio ln_cons_exp_all ///
         urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year i.region_num, q(`q')
    
    * Save regression results for each quantile
    estimates store q`q'
}

esttab q10 q20 q30 q40 q50 q60 q70 q80 q90 using "found_num.rtf", ///
   b(3) se(3) star(* 0.1 ** 0.05 *** 0.01) ///
   nogaps compress ///
   stats(year_fe region_fe N r2_p, ///
       labels("Year FE" "Region FE" "Observations" "Pseudo R2") ///
       fmt(%9.0g %9.0g %9.0g %9.3f)) ///
   mtitles("Q10" "Q20"  "Q30"  "Q40" "Q50"  "Q60"  "Q70" "Q80" "Q90") ///
   keep(treatment ln_direct_donation total_comm_num ///
        labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve) ///
   replace title("Quantile Regression Results-found_num") ///
   rtf


foreach q in 10 20 30 40 50 60 70 80 90 {
    qreg2 social_org_num treatment ln_direct_donation  ///
         total_comm_num  labor_gdp_ratio ln_cons_exp_all ///
         urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year i.region_num, q(`q')
    
    * Save regression results for each quantile
    estimates store q`q'
}

esttab q10 q20 q30 q40 q50 q60 q70 q80 q90 using "social_org_num.rtf", ///
   b(3) se(3) star(* 0.1 ** 0.05 *** 0.01) ///
   nogaps compress ///
   stats(year_fe region_fe N r2_p, ///
       labels("Year FE" "Region FE" "Observations" "Pseudo R2") ///
       fmt(%9.0g %9.0g %9.0g %9.3f)) ///
   mtitles("Q10" "Q20"  "Q30"  "Q40" "Q50"  "Q60"  "Q70" "Q80" "Q90") ///
   keep(treatment ln_direct_donation total_comm_num ///
        labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve) ///
   replace title("Quantile Regression Results-social_org_num") ///
   rtf


foreach q in 10 20 30 40 50 60 70 80 90 {
    qreg2 priv_npe_num treatment ln_direct_donation  ///
         total_comm_num  labor_gdp_ratio ln_cons_exp_all ///
         urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year i.region_num, q(`q')
    
    * Save regression results for each quantile
    estimates store q`q'
}

esttab q10 q20 q30 q40 q50 q60 q70 q80 q90 using "priv_npe_num.rtf", ///
   b(3) se(3) star(* 0.1 ** 0.05 *** 0.01) ///
   nogaps compress ///
   stats(year_fe region_fe N r2_p, ///
       labels("Year FE" "Region FE" "Observations" "Pseudo R2") ///
       fmt(%9.0g %9.0g %9.0g %9.3f)) ///
   mtitles("Q10" "Q20"  "Q30"  "Q40" "Q50"  "Q60"  "Q70" "Q80" "Q90") ///
   keep(treatment ln_direct_donation total_comm_num ///
        labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve) ///
   replace title("Quantile Regression Results-priv_npe_num") ///
   rtf


* Conduct quantile regression analysis
foreach q in 10 20 30 40 50 60 70 80 90 {
    qreg2 total_social_org treatment ln_direct_donation  ///
         total_comm_num  labor_gdp_ratio ln_cons_exp_all ///
         urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year i.region_num, q(`q')
    
    * Save regression results for each quantile
    estimates store q`q'
}

esttab q10 q20 q30 q40 q50 q60 q70 q80 q90 using "total_social_org.rtf", ///
   b(3) se(3) star(* 0.1 ** 0.05 *** 0.01) ///
   nogaps compress ///
   stats(year_fe region_fe N r2_p, ///
       labels("Year FE" "Region FE" "Observations" "Pseudo R2") ///
       fmt(%9.0g %9.0g %9.0g %9.3f)) ///
   mtitles("Q10" "Q20"  "Q30"  "Q40" "Q50"  "Q60"  "Q70" "Q80" "Q90") ///
   keep(treatment ln_direct_donation total_comm_num ///
        labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve) ///
   replace title("Quantile Regression Results-total_social_org") ///
   rtf
  
 
 // Spatial Econometrics: Convert GIS map data to dta format - Province_shp.dta and Province.dta
preserve
spshape2dta 省, replace
use 省.dta, clear 
rename ʡ____ fips
save "省.dta", replace
restore
 
 * Assume the current dataset has a variable 'region' corresponding to region names
 * First generate a new variable 'fips'
gen fips = .

* Match and replace one by one according to the mapping table:
replace fips = 110000 if region == "北京"
replace fips = 120000 if region == "天津"
replace fips = 130000 if region == "河北"
replace fips = 140000 if region == "山西"
replace fips = 150000 if region == "内蒙古"
replace fips = 210000 if region == "辽宁"
replace fips = 220000 if region == "吉林"
replace fips = 230000 if region == "黑龙江"
replace fips = 310000 if region == "上海"
replace fips = 320000 if region == "江苏"
replace fips = 330000 if region == "浙江"
replace fips = 340000 if region == "安徽"
replace fips = 350000 if region == "福建"
replace fips = 360000 if region == "江西"
replace fips = 370000 if region == "山东"
replace fips = 410000 if region == "河南"
replace fips = 420000 if region == "湖北"
replace fips = 430000 if region == "湖南"
replace fips = 440000 if region == "广东"
replace fips = 450000 if region == "广西"
replace fips = 460000 if region == "海南"
replace fips = 500000 if region == "重庆"
replace fips = 510000 if region == "四川"
replace fips = 520000 if region == "贵州"
replace fips = 530000 if region == "云南"
replace fips = 540000 if region == "西藏"
replace fips = 610000 if region == "陕西"
replace fips = 620000 if region == "甘肃"
replace fips = 630000 if region == "青海"
replace fips = 640000 if region == "宁夏"
replace fips = 650000 if region == "新疆"

merge m:1 fips using "省.dta"
drop if _merge==2
drop _merge
duplicates drop fips year , force

xtset _ID year
spset

// Create spatial weight matrix

spmatrix create contiguity M3 if year== 2016, norm(row) second(0.5) replace 
// Create row-standardized spatial weight matrix
// norm(row) indicates row standardization
// second(0.5) sets the weight for second-order neighbors to 0.5

* Part 2: Export matrix to text file
spmatrix export M3 using contig2.txt ,replace    // Export row-standardized matrix


* Part 3: Read and process each weight matrix file

* Process contig2.txt
preserve
clear
insheet using"contig2.txt", delimiter(" ")
drop if _n==1  
drop v1
rename v# v#, renumber sort  
save contig2.dta, replace
restore

* Part 4: Create standardized spatial weight matrix
spatwmat using contig2 , name(W3) standardize    // Row-standardized first and second-order contiguity spatial weight matrix W3

// Spatial Autocorrelation Test
// Row-standardized first and second-order contiguity spatial weight matrix W3
forvalues y = 2006/2016 {
    preserve  // Save current data state
    keep if year == `y'
    display _newline "===== Results for Year `y' ====="
    spatgsa total_social_org  , weights(W3) moran
    restore  // Restore data state
}


* ln_direct_donation  c_so_don_cash_yi c_social_don_total_yi
* Model 1: Base model  found_num social_org_num priv_npe_num
xsmle total_social_org treatment ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(treatment) model(sdm) nolog

xsmle found_num treatment ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(treatment) model(sdm) nolog

xsmle social_org_num treatment ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(treatment) model(sdm) nolog

xsmle priv_npe_num treatment ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(treatment) model(sdm) nolog


* Model 2: Treatment × social_org_num - Analysis of interaction between treatment effect and direct donations

preserve
keep if inrange(year, 2010, 2015)
xsmle total_social_org c.c_treatment##c.c_so_don_cash_yi ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(c.c_treatment##c.c_so_don_cash_yi) model(sdm) nolog
restore

preserve
keep if inrange(year, 2010, 2015)
xsmle found_num c.c_treatment##c.c_so_don_cash_yi ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(c.c_treatment##c.c_so_don_cash_yi) model(sdm) nolog
restore

preserve
keep if inrange(year, 2010, 2015)
xsmle social_org_num c.c_treatment##c.c_so_don_cash_yi ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(c.c_treatment##c.c_so_don_cash_yi) model(sdm) nolog
restore

preserve
keep if inrange(year, 2010, 2015)
xsmle priv_npe_num c.c_treatment##c.c_so_don_cash_yi ln_direct_donation  total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(c.c_treatment##c.c_so_don_cash_yi) model(sdm) nolog
restore

* Model 3: Treatment × Government management expenditure interaction 
xsmle total_social_org c.c_treatment##c.c_ln_ngo ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(c.c_treatment##c.c_ln_ngo) model(sdm) nolog 

xsmle found_num c.c_treatment##c.c_ln_ngo ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(c.c_treatment##c.c_ln_ngo) model(sdm) nolog 

xsmle social_org_num c.c_treatment##c.c_ln_ngo ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(c.c_treatment##c.c_ln_ngo) model(sdm) nolog 

xsmle priv_npe_num c.c_treatment##c.c_ln_ngo ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(c.c_treatment##c.c_ln_ngo) model(sdm) nolog 

preserve
keep if inrange(year, 2010, 2015)
xsmle total_social_org c.c_treatment##c.c_social_don_total_yi ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(c.c_treatment##c.c_social_don_total_yi) model(sdm) nolog
restore

preserve
keep if inrange(year, 2010, 2015)
xsmle found_num c.c_treatment##c.c_social_don_total_yi ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(c.c_treatment##c.c_social_don_total_yi) model(sdm) nolog
restore

preserve
keep if inrange(year, 2010, 2015)
xsmle social_org_num c.c_treatment##c.c_social_don_total_yi ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(c.c_treatment##c.c_social_don_total_yi) model(sdm) nolog
restore

preserve
keep if inrange(year, 2010, 2015)
xsmle priv_npe_num c.c_treatment##c.c_social_don_total_yi ln_direct_donation total_comm_num labor_gdp_ratio ln_cons_exp_all urban_rate gdp_growth private_emp_ratio ln_pop ln_srve i.year, fe wmat(W3) durbin(c.c_treatment##c.c_social_don_total_yi) model(sdm) nolog
restore